The mass of the 6-quark from lattice NRQCD and lattice perturbation theory 
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We present a determination of the 6-quark mass accurate through 0{al) in perturbation theory 
and including partial contributions at 0{as)- Nonperturbative input comes from the calculation 
of the T and energies in lattice QCD including the effect of u, d and s sea quarks. We use an 
improved NRQCD action for the b-quark. This is combined with the heavy quark energy shift in 
NRQCD determined using a mixed approach of high-/? simulation and automated lattice pertur- 
bation theory. Comparison with experiment enables the quark mass to be extracted: in the MS 
scheme we find ffibijnb) ~ 4.166(43) GeV. 

PACS numbers: 12.38.Bx, 12.38.Gc 



I. INTRODUCTION 

The accurate determination of quark masses is an im- 
portant component of high-precision tests of the Stan- 
dard Model. Because quarks cannot be isolated experi- 
mentally, the mass must be defined carefully and its ex- 
traction from quantities that are accessible to experiment 
must be well controlled from the theory side. The 6-quark 
mass is particularly important: its uncertainty feeds into 
errors in tests of the Standard Model in B physics as well 
as into the cross-section for the Higgs decay, H — > bb. 

The most accurate results to date for the 6-quark mass 
come from comparison of the experimental cross section 
for e+e~ to hadrons in the bottomonium region with 
high-order (af ) continuum QCD perturbation theory [T]- 
|3]. Errors of 0.5% are possible. A similar method has 
now been applied to lattice QCD results [H [S], using 
pseudoscalar correlators made from heavy quarks instead 
of the experimental cross-section. For these calculations, 
the experimental input is the value of the meson mass (in 
this case the r](,) used to tune the lattice 6-quark mass. 
Again a 0.5% error is achieved and good agreement is 
seen with the continuum results. 

It is important to test these determinations against 
a different method of obtaining the 6-quark mass which 
has completely uncorrelated systematic errors. This is 
the aim of this paper. We use a direct determination 
from full lattice QCD calculations of the binding energy 
of both T and Bg mesons. Since we use a nonrelativistic 
effective theory for the 6-quark (NRQCD) [Bl[7] this needs 
a calculation of the heavy quark energy shift. We do this 
in lattice QCD perturbation theory through two-loops 
(with partial three- loop contributions), significantly im- 
proving on earlier determinations that used one-loop cal- 
culations [5]. We have also implemented a one- loop im- 
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proved NRQCD action to reduce systematic errors. 



Calculating higher order loop corrections in lattice per- 
turbation theory for heavy quarks in NRQCD grows ever 
more difficult with each order owing to the increasing 
number of diagrams and the complicated vertex struc- 
ture. Various authors [QlflT] have suggested an approach 
in which the heavy quark propagator is measured in the 
weak coupling regime and the renormalization param- 
eters are fitted to a polynomial in a^, thus obtaining 
the radiative corrections beyond one loop. This method 
is certainly practical for obtaining the quenched contri- 
butions to renormalization parameters since quenched 
gauge configurations are relatively cheap to generate. At 
two loop order there are relatively few remaining dia- 
grams with sea quark loops and these can be feasibly 
computed using automated lattice perturbation theory. 
In contrast, there are many two loop diagrams contain- 
ing only gluon propagators that pose a challenging task 
for direct evaluation with automated lattice perturbation 
theory. We therefore employ a mixed approach to the 
determination of the two-loop heavy quark energy shift 
combining quenched high-/3 calculations with automated 
lattice perturbation theory for the sea quark pieces. 



In Section |TT] we discuss how we extract the 6-quark 
mass from simulations of lattice NRQCD. Section [ill A| 
describes the automated lattice perturbation theory com- 
putation of the fermionic contributions to the two loop 
energy shift. We present our implementation of the high- 
/3 method in Section |IIIB| including the concomitant fi- 
nite volume perturbation theory in Appendix [X] The de- 
tails of the standard non-perturbative part of the cal- 
culation are given in Section |IV[ Finally we detail the 
extraction of the MS mass in Section |v] and present our 
conclusions in Section IVlIl 
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II. EXTRACTING THE 6-QUARK MASS 



B. Matching the pole mass to the MS mass 



Quark confinement ensures that quark masses are not 
physically measurable quantities, so the notion of quark 
mass is a theoretical construction. A wide range of 
quark mass definitions exist, often tailored to exploit the 
physics of a particular process. One common choice of 
quark mass is the pole mass, defined as the pole in the 
renormalised heavy quark propagator. The pole mass, 
however, is a purely perturbative concept and suffers 
from infrared ambiguities known as renormalons ^12} I13j . 
A better mass is the running mass in the MS scheme, 
which is free of renormalon ambiguities by construction, 
and is the usual choice for quoting the quark masses. 
Lattice calculations use the renormalon-free bare lattice 
mass which must then be matched to MS to enable 
meaningful comparison. We match bare lattice quan- 
tities to the MS mass using the pole mass as an inter- 
mediate step. Any renormalon ambiguities cancel in the 
full matching procedure between the lattice quantities 
and the MS mass, as we argue below. For an explicit 
demonstration, see [Ml. 



A. Extracting the pole mass 

We determine the heavy quark pole mass, M^oie, by 
relating it to the experimental T mass M™^*. The mass 
of a heavy meson is given by twice the pole quark mass 
plus the binding energy. In an effective theory such as 
NRQCD, physics above the scale of the 6-quark mass is 
removed and the zero of energy for the heavy quark is 
shifted by Eq, leading to the relation |15]: 



2M 



pole 



(1) 



A/po/e = M^^^P* - a-\aE^^^^ - aE^) 



Here i^sim is the energy of the T meson at zero momen- 
tum, extracted from lattice NRQCD data at lattice spac- 
ing a. The quantity (-Esim ~ Sii'o) corresponds to the 
"binding energy" of the meson in NRQCD and we must 
determine Eq perturbatively in order to find Mpoie- With 
our NRQCD action, we can also calculate the pole mass 
using the meson 

(2) 

We use this as a check for systematic errors which could 
be quite different in heavy-heavy and heavy-light sys- 
tems. 

In principle one could extract the quark mass by di- 
rectly matching the pole mass to the bare lattice NRQCD 
mass in physical units, toq, via the heavy quark mass 
renormalisation, , 

Mpoie = Zmo{ami3)mo. (3) 

We found, however, that extracting a sufficiently pre- 
cise quenched two loop mass renormalisation from high-/? 
simulations was not possible with the statistics available. 
In this paper, we therefore discuss only the energy shift 
method. 



The mass renormalisation relating the pole mass to the 
MS mass, mf,, evaluated at some scale is given by 



r71fc(/i) = Zj^l{^i)Mpole, 



(4) 



and has been calculated to three- loops in [IB] . 

Although the pole mass is plagued by renormalon am- 
biguities, these ambiguities cancel when lattice quantities 
are related to the MS mass. This can be seen by equat- 
ing Eqns ([T]) and ^ and rearranging them to obtain 

-^niQ- Eq) = Af^'^P* 



2(Z„ 



(5) 



The two quantities on the right hand side of the equa- 
tion are renormalon ambiguity free: M^'^p* is a physical 
quantity and £'sim is determined nonperturbativcly from 
lattice simulations. Any renormalon ambiguities in the 
two power series, Z^o and Eg, on the left-hand side of 
the equation must therefore cancel at every order in as- 
This renormalon cancellation is also evident in the direct 
matching of the bare lattice mass to the MS mass, 

mtifi) = Z,„„(amo)Zj^/(/i)mo, (6) 

as both ffib and mo are renormalon-free. 

We combine Eqns ([T]) and ^ to relate lattice quanti- 
ties to the MS mass 

fribifi) = ^Z]^/{p) [M^'^P* - a-\aEsi^ - 2aEo)] , (7) 

and similarly for the Bg meson 

mhip) - [Af]j"/* - a-\aE,;„,^B, - aEo)] . (8) 

These relations will be used to extract mf,(mb) once we 
have calculated Eq and i?sim, which we describe in detail 
in the next sections. 



C. NRQCD, gluon and light quark actions 

We now describe the heavy quark, gluon and light 
quark actions used in our calculation. We use the 
Symanzik improved 0{v*) NRQCD action, given in 
[H |T7], which has already been successfully used by 
HPQCD in a number of heavy quark physics calculations, 
see e.g. [51 [TT^ET] . The Hamiltonian is given by 



aH = aHo + aSH; 



i5H 



2amo 

(A(2))2 



(9) 
(10) 



-ci 



-C3 



-C6 



8(amo)'^ 

9 ^ 
8(amo)2 

2amo 
a(A(2))2 
16n(amo)2 



C2 



«3 



8(amo)2 
V X E - E X V 



V-E - E-V 



C5 



g^AW 
24amo 



(11) 
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TABLE I: Values of the one loop corrections in the 
Ci = 1.0 + QsC^^' at two bare masses, and the scale at 
each coefficient is evaluated. 



where 



series 
which 



CoefScient 


amo = 2.5 


amo = 1.72 


9* 


Cl 


0.95 


0.766 


1.8/a 


C4 


0.78 


0.691 


Tr/a 


C5 


0.41 


0.392 


1.4/a 


C6 


0.95 


0.766 


1.8/a 



A(2),V and AW are covariant lattice derivatives, E 
and B are improved chromo-electric and magnetic field 
strengths, n is a stability parameter that will be de- 
scribed below and amo is the bare &-quark mass in lattice 
units. The are the Wilson coefficients of the effective 
theory and the terms are normalised such that they have 

the expansion q = l+asc'f'^ +0{a'^). All gauge fields are 
tadpole improved with the fourth root of the plaquette 
uo,p- 

The one loop corrections c^'^ are described in [17] and 
we include these for ci, C4, C5, Cg in the high-/? simulation 
and the nonperturbative determination of -Esim- The c^^"^ 
are a function of the effective theory cutoff, in this case 
the bare quark mass amo, but the total coefficient will 
also depend on the scale for as- We estimate the ap- 
propriate scale for several of the coefficients using the 
BLM procdure [22| which gives q* = 1.8/a for ci,cq and 
q* = 1.4/a for C5. For C4 we take q* = n/a. The values 
of the one loop corrections for two bare masses relevant 
to this calculation are given in table [l] We use as in the 
l^-scheme. 

The 6-quark propagators are generated by time evolu- 
tion using the equation 



(3pi = 



Ppg — 



10 

7^' 



2 (1 + 0.4805a,). 



0.03325a.,. 



(14) 
(15) 
(16) 



Mo.p is the tadpole improvement factor coming from the 
fourth-root of the plaquette. The same action is used 
for the MILC gauge configurations used in the non- 
perturbative determination of i?sim and for the high-/? 
simulations. The action in the high-/? simulations in- 
cludes an additional factor coming from the use of twisted 
boundary conditions, see section |IIIB[ The value of a, 
used in the improvement coefhcicnts is given by the for- 
mula used by the MILC collaboration [25j 



as = 1.3036 log(uo,p(/3)). 



(17) 



Here we use the quenched values of uq_p{/3) determined 
from our high-/? configurations. The MILC configura- 
tions used in our nonperturbative analysis include sea 
quarks and so have additional 0{nfa^) contributions. 
However, these only affect Eq at 0(n/a^) and so appear 
in terms we have not calculated anyway. These terms 
are part of our error budget. We give more details of the 
generation of high-/? configurations in Appendix [B| 

Light sea quarks are included with the ASQtad im- 
proved staggered action [55] in both the nj — 2 + 1 MILC 
gauge configurations used to determine Egim [23 HZ] and 
in the automated perturbation theory for Eq. 



III. PERTURBATIVE DETERMINATION OF 
THE HEAVY QUARK ENERGY SHIFT 



Gix,t + 1) ^ 1 



aSH 
2n 



2n 
a6H 



G{x,t) (12) 



Here we first describe the calculation of the one-loop 
contribution and the two-loop fermionic contribution to 
Ep. The high-/? method used to compute the gluonic 



two-loop contribution is described in the section HI B 



for some initial condition G{x^Q). The parameter n is 
included for numerical stability and is set to 4, which is 
sufficient for all quark masses used here. Once it is high 
enough, results do not depend on the value of n [5]. 

The gluon action is a Symanzik improved Liischer- 
Weisz action [23l [24] 

Slw[U] = /?p,^-^ReTr(l-{/pO 

X 

+ /3rt^^RcTr(l-C/rt) 

+ /3psE]^I^cTr(l-[/p,), (13) 



A. Automated lattice perturbation theory 

We calculate the one loop gluonic and the two loop 
sea quark contributions to the heavy quark renormaliza- 
tion constants using the automated lattice perturbation 
theory routines HiPPy and HPSRC ^28,^. These rou- 
tines have now been widely used and extensively tested 
in a variety of perturbative calculations, for example in 
[ID, 17, 30. .35]. 

Evaluating the relevant Feynman integrals with 
HiPPy and HPsRC is a two-stage process: firstly the 
python routine HiPPy generates Feynman rules encoded 
in "vertex files" . These vertex files are then read in by the 
HPSRC code, a collection of FORTRAN modules that 



4 



reconstruct the diagrams and evaluate the corresponding 
integrals numerically, using the VEGAS algorithm (36] , 
All derivatives of the self energy are implemented ana- 
lytically using the derived taylor type, defined as part 
of the TaylUR package ^7.. Our computations of these 
diagrams were performed on the Darwin cluster at the 
Cambridge High Performance Computing Service with 
routines adapted for parallel computers using Message 
Passing Interface (MPI). 

There are several advantages associated with using 
automated lattice perturbation theory, and the HiPPy 
/HPSRC routines in particular. First, automation re- 
moves the need to manipulate complicated expressions 
by hand. Secondly, the modular nature of the HiPPy 
and HPSRC routines greatly simplifies the use of differ- 
ent actions. Once Feynman diagrams are encoded in an 
HPsRC routine, the same calculation can be easily re- 
peated with different quark and gluon actions by simply 
changing the input vertex files. This allows one to rel- 
atively easily reproduce previously published results for 
different actions, which serves as a nontrivial check of the 
routines. 

Furthermore, the modules in HPSRC can be reused. 
We took advantage of this for the two loop calculations 
presented in this paper: the same fermionic insertions in 
the gluon propagator appear in the two loop diagrams 
for both the heavy quark energy shift and the tadpole 
improvement factor, uq. 

We wrote two "skeleton" one loop HPSRC routines: 
one to calculate the one loop energy shift and one for 
the one loop tadpole improvement factor. Reproducing 
previously published results, such as those in [SF and 
[55] respectively, confirmed that these one loop routines 
were correct. The corresponding two loop diagrams (see 
Figure [T]) are simply the one loop skeleton diagrams with 
the "bare" gluon propagator replaced by the "dressed" 
gluon propagator that includes the fermion insertions; 
these insertions were calculated in a separate routine 
gluon_sigma. This routine was debugged by confirm- 
ing that the appropriate Ward identity was satisfied by 
the dressed gluon propagator. 

At two loops there are four diagrams with internal 
fermions that contribute to the energy shift. We il- 
lustrate these contributions in Figure [T] Double lines 
are heavy quark propagators coming from the improved 
NRQCD action, single lines are ASQtad sea quark prop- 
agators and curly lines are from the Symanzik improved 
gluon action. The radiative corrections to the NRQCD 
and ASQtad actions described in section [II C| are not in- 
cluded in the perturbative calculation as these only affect 
Eq at higher order in a^. 

We calculated the heavy quark energy shift at two dif- 
ferent heavy quark masses discussed in section |IV[ At 
each heavy quark mass we use nine different light quark 
masses and extrapolate to zero light quark mass. We 
tabulate our extrapolated results in Table [V[ where they 

— (2) 

appear as the n /-dependent contribution to Eq ' . 

The energy shift is infrared finite, but we introduce a 




FIG. 1: Fermionic contributions to i?o, calculated using au- 
tomated lattice perturbation theory. Double lines indicate 
heavy quarks, curly lines are gluons and single lines represent 
light sea quarks. 

gluon mass as an intermediate regulator to ensure con- 
vergence for the numerical integration. We confirmed 
that the results are independent of the gluon mass for 
sufficiently small gluon mass, which in this case was ap- 
proximately c?)? < 10~^. 

Wc will also need the sea quark contribution to the 
tadpole improvement factor uq since the high-/? simula- 
tion includes only the gluonic piece. We calculate this 
using the automated perturbation theory. The perturba- 
tive expansion for the tadpole factor is written as 

u„ = l-ul''>aL-ul^^al+0{al). (18) 

The two loop expansion for the plaquette tadpole is given 
by Mason ^\ and we explicitly computed the one-loop 
coefficient and the two- loop n f coefficient which we quote 
here and which both agree with Mason. The result is 

uo,p = 1 - 0.76708(2)aL 

" (1.7723 - 0.069715(7)?i/)a^ + 0{al). (19) 

We require only the coefficent of n/a|. For completeness 
we also computed the two-loop nj contribution to the 
Landau tadpole. The quenched two-loop Landau tadpole 
was computed by Nobes et al. [39] and together with our 
result the Landau tadpole is 

wo,L = 1 - 0.7501(l)aL 

- (2.06(1) - 0.0727(l)n/)a| + 0{al). (20) 

B. The high-/? method 

The high-/? method allows us to compute the gluonic 
contributions to the quark propagator by generating an 
ensemble of quenched lattice gauge configurations at very 
weak coupling and calculating the dressed 6-quark prop- 
agator. The energy of the propagator can then be de- 
scribed very well by a power series in the QCD coupling, 
which we fit to the Monte-Carlo data to extract the rel- 
evant two-loop and higher contributions to Eq. 

It is important in high-/3 studies to eliminate non- 
perturbative contributions which are due to the tun- 
nelling of fields and their associated Polyakov lines, or 
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torelons, between vacua associated with toron gauge 
configurations [41]. Such tunneUing is suppressed using 
twisted boundary conditions [I^l - Hi] for which there is no 
zero mode for the non-abeUan gauge field. The Polyakov 
Hue that traverses aU the directions with twisted bound- 
ary conditions has a non-zero expectation value for a 
given configuration. This expectation value is complex 
and if no tunnelling has occured it is proportional to 
an element of Z3. We verify that this is the case for 
the configurations we use. As is shown later in this sec- 
tion, see the discussion leading to Eqs. (43) (44), twisted 



boundary conditions also considerably reduce finite-size, 
L-dependent effects which significantly aids the fitting 
process. 

We carry out the high-/? simulation on finite size lat- 
tices of volume x T, with typically T = 3L, for 
a range of values for (3 and L. Here L is the spatial 
extent and T the temporal extent of the lattice. We 
use L values from 3 to 10 inclusive and /3p/ values of 
12,15,16,20,24,27,32,38,46,54,62,70,80,92, and 120. We 
then perform a simultaneous fit in as and L to deduce 
the L — > 00 limit for the expansion of measured quanti- 
ties as a power series in as- 

We denote the gauge fields by Uf^{x), and on a lattice 
with sites in the /x direction they satisfy the boundary 
condition 

where the twist matrices are defined by 



(21) 



z — exp {27:i/Nc) , e (0,...,iV,-l) .(22) 

Here n^i/ is antisymmetric and its values must be chosen 
so that t^ivapn^unap — 0\n^- This choice ensures config- 
urations have zero topological charge. For Nc = 3 we 
apply a non-trivial twist in the spatial directions, which 
we label 1, 2 and 3, with ni2 = nis = 7123 = 1 and 
= 0. 

With twisted boundary conditions, the fermion fields, 
ip, are Nc x Nc colour-times-smcU matrices. "Smell" is 
a new quantum number that allows twisted boundary 
conditions to be applied to fermion fields; colour labels 
the rows and smell the columns. Then, as for the gauge 
fields, 



(23) 



Under a gauge transformation given by the SU{Nc) 
field g{x) the quantum fields transform as 



Upix) g{x)Uf,{x)g\x + Ef,), 



(24) 



where g{x + Lye^) — ri^g{x)ril. We define the auxiliary 
gauge fields 



L,. 



(25) 



Then under a gauge transformation U^{x) transforms as 
in Eq. (24) but now with g{x) regarded as periodic: g{x+ 



L^,e^) = g{x). 

The gauge action is of the form 



S{U)^P J2 cpfp{x)P(U,x) 



(26) 



where A is the set of all lattice sites; P(/7, x) is the trace 
over a general Wilson loop; cp is a numerical coefficient 
and fp{x) £ Zjv^ is a phase factor defined by 



!p{x) = n (^"^'')" 



(27) 



p<v 



Here uj^^{P,x) is the winding number of the Wilson 
loop projected onto the (fj,, v) plane about the point 
x^ = x^ = (L + 1/2). An explicit representation for 
the twist matrices fi^ is not needed to compute fp{x). 
When fermions are included, however, the implementa- 
tion of twisted boundary conditions for general Wilson 
lines does require a representation for the Sl^ to be cho- 
sen. 

One method for implementing the boundary conditions 
extends the lattice by tiling with twisted periodic transla- 
tions of the original configuration, effectively surround- 
ing the lattice with a halo of links. This method has 
major disadvantages: it is difficult to parallelize because 
the physical sites are a subset of the tiled lattice array; 
it requires more storage; and in improved NRQCD the 
Wilson lines can extend far into the tiled region, which 
means that the extent of the halo needs to be significant. 
Rather than extending the lattice we write the action 
in terms of the auxiliary gauge fields, U^{x). Then one 
can show that all Wilson lines can be constructed using 
the auxiliary gauge fields with periodic boundary condi- 
tions multiplied on the right by an SU{Nc) matrix. This 
SU{Nc) matrix, which we denote R{V), is constructed 
from a product of the twist matrices, fl^, and is deter- 
mined by the ordered and signed sequence in which the 
line crosses the boundaries. We now discuss this con- 
struction in more detail. 

A general path V{x,y;s) starting at site a; on a lat- 
tice in dimension D is defined by an ordered list s = 
[sq, Si, . . . Si_i] of signed integers, s^, 1 < \si\ < D, which 
denote the steps along the path. The j-th point on the 
path is Zj where 



zq = X, Zj+i ^ Zj +es. < j <l , 



(28) 



with endpoint defined hy y — zi. We define the ordered 
product of links along the path Vix, y; s) to be 



C{U;V) 



tY[UsM) 



i=0 



where, for ^ e {1, 2, . . . , D}, 

U^,{x)^Ul{x-e,), 



(29) 



(30) 
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and the U fields satisfy the periodic boundary condition 



(31) 



The ordering operator T means that matrices in the 
product are ordered from left to right with increasing 
index, i. The Wilson line L{U;'P) associated with the 
path V{x, y; s) is then 



L{U\V)= C{U;V)R{V) . 



(32) 



To implement the twisted boundary conditions without 
using a lattice halo we define the SU{Nc) matrix as fol- 
lows. A list [cq{x), Ci{x), . . . , Cp-i{x)] is associated with 
the Wilson line starting at x, where the Cj are signed 
integers 1 < \cj\ < D. The line crosses a boundary of 
the hypercube p times. On the j-th crossing it crosses a 
boundary in a direction parallel to the fj,j axis in the pos- 
itive (negative) direction. We define the corresponding Cj 
to be Cj — ~fij{iij). R{P) is then given by 



RiV) 




(33) 



with the convention fi_^ = il^ and where T is the index- 
ordering operator defined above. L{U;V) is then the 
parallel transporter from the endpoint y back to the start 
point X. By expressing the Wilson line in terms of the U 
fields the boundary conditions are implemented simply 
by right-multiplication by RiV). A similar result holds 
for the evolution of the NRQCD Green function as we will 
describe below. With these conventions, a Wilson loop 
W{x, s), located at x and defined by the path 'P{x, x; s), 
is given by 



W{x,s) 



1 

n: 



TT(LiU;V{x,x;s)) 



(34) 



The basis states for the fermion field ip are the 
independent x N^., colour-times-smell, real matrices. 
Twisted boundary conditions admit fractional momenta 
on the lattice and for twisted boundary conditions in the 
1, 2, 3 directions and periodic boundary conditions in the 
fourth direction, the allowed momenta are of the form 



P 
k 



(ni, 712,^3,0) -I- fc. 



27r 

ttIi Trl2 ttIs ttI^ 

IT' IT' IT' 'Y 



(35) 



where the Z^, for r = 1,2,3, are integers with — L/2 < 
Ir < L/2 and -T/2 < ^4 < T/2 {L and T assumed 
even). The possible entries in the integer vector n = 
(ni, n2, na, n4) depend on the number of directions in 
which the boundary condition is twisted. In our case we 
have < ni,n2 < Nc, ns = — (ni -|- 'T-2)|Ar^ and n4 = 0. 



In NRQCD the source on the initial time slice for a 
Green function with momentum p is 

1 



X{p,x) 



ipx 



(ni+n2){ni+n2- 



1) n^""^ 




(36) 

We need an explicit representation for the O^, where /i = 
1, 2, 3, and for iVc = 3 we choose 

/o 1 o\ 

02= 1 n\nl . (37) 
\l o) 

In the case of purely periodic boundary conditions we 
can take the source for the NRQCD Green function to 
be 1 • e'P'^, where 1 is the Nc x Nc unit matrix. This 
evolves all quark colour states in one go. The analogous 
approach for quarks labelled by colour times smell is not 
convenient and so we evolve a source appropriately cho- 
sen from the basis of Nc x Nc matrices described above; 
colour and smell singlet states, if needed, must then be 
constructed explicitly. The Green function G{x,p, t) sat- 
isfies the usual twisted boundary conditions 

G{x + Le^,p,t)^Vl^G{x,p,t)^l . (38) 

The NRQCD evolution for G{x,p,t) is given by the 
full NRQCD action and takes the form 



Gix,p,t + l) = }2K{x,y,t)G{y,p,t) 
y 

The kernel K is given by 



(39) 



Kix,y,t) = 1- 



SH 



1 - 



Ho 



2n 

Ho 
2n 



SH 



(40) 



with iJo, SR defined in section [IIC[ 

We implement the operators in K using a python pre- 
processing package that defines each operator in Ho and 
SH as a list of Wilson paths. The Wilson paths are each 
defined by a list s with a complex amplitude; these op- 
erator definitions are read in at run time. We apply the 
action of each operator on G{y,p,t) with a standard 
function that: first constructs the parallel transporter 
L{U ; 'P{x, y s)) for each path weighted by the associated 
amplitude where x = {x,t -t- l),j/ = {y,t), then per- 
forms the parallel transport of G from the t-ih to the 
[t + l)-th time slice, and finally accumulates the results 
in G{x,p,t + 1). We solve the problem of implement- 
ing the twisted boundary conditions in carrying out this 
calculation by using U fields. The net result is that the 
evolution equation can be written as 



G{x,p,t + l) = J2bm (^/:(C/;P™) 

m \ y 



y.G{y,p,t)R{V^)) , (41) 
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where Vm = 'P{x, y, Sm) and the sum over m runs over 
all lists, Sm, that define the kernel K{x,y,t), with bm 
the amplitude of the m-th line. The matrix R{Vm) im- 
plements the twisted boundary conditions and is simple 
to compute for each Vm- Because R right- multiplies the 
Green's function and time evolution is a left-multiplying 
operation, we can perform the time evolution for a given 
m using periodic boundary conditions for the U fields 
and then independently right-multiply by the associated 
R matrix. This method removes the need for any halo 
of gauge fields and the whole calculation can be easily 
parallelized. 

Furthermore, twisted boundary conditions reduce 
finite-size effects in colour singlet observables. To illus- 
trate this result, we can consider the example of the cor- 
relator for a meson at rest, which is given by 



M{t) 



E 

ya 



Tr[G„(y,0;t)Gl(y,0;0], 



(42) 



where a labels the basis matrix used for the source of the 
quark propagator located at the origin; all irrelevant spin 
degrees of freedom have been suppressed. The correlator 
M{t) is the sum of weighted Wilson loops consisting of 
a Wilson line Li, connecting x = (0,0) to y = {y,t), 
followed by L2 connecting y back to x and defined by the 
paths Vi = 'P{x, y, Si) and V2 = 'P[y^ x, S2), respectively. 
Then M{t) is of the form 

Af(t)= J2 fiVi,V2)Tv[C,XaRiR2xiC2], (43) 

'Pl,'P2,y,a 

where /(7^i,7^2) is the amplitude associated with the 
loop, {Vi +V2), A = CiU;V^) and Ri = R(r^), for 
i — 1,2. Irrespective of the details of £1 and £2, the 
term sandwiched in the middle is 

;Xai?ii?2xl =Tr[i?ii?2]l. (44) 



Since R1R2 is a product of the 17 matrices and their con- 
jugates, the trace in the above formula vanishes unless 
R1R2 = 1. Thus, for a non-zero contribution, the Wil- 
son loop composed of Vi and 7^2 must loop around the 
spatial torus a multiple of Nc times in such a way that 
R1R2 = 1. This reduces finite size effects as the effec- 
tive size of the lattice is now of order N^. times its spatial 
extent. 



1. Perturbative fitting of Eq 

We obtain the quark propagator by averaging G{x, p, t) 
over the ensemble of high-/? configurations. Because 
G{x,p,t) is not gauge invariant we fix the configura- 
tions to Coulomb gauge. We then define the Coulomb 
ensemble-averaged quark propagator by 



Here we write G(p, t, /?, L) as a function of L to indicate 
explicitly that there are finite size effects, which must be 
accounted for to extract the desired L — > 00 result. 

In order to extract the two-loop and three-loop co- 
efficients in the perturbation expansion for Eq using the 
high-/? method it is necessary to carry out a simultaneous 
two parameter fit in and L. The fit is a power series 
in as and in 1 /L and we measure the L — >■ 00 coefficient 
of the a", for n — 2,3, terms. Because the signal for the 
two- loop, Qfg, term is small compared with the one- loop 
contribution the accuracy of the fit is greatly improved 
by calculating the one-loop coefficient analytically and 
so determining the coefficient of as in the fit. However, 
Feynman perturbation theory on the lattice gives the re- 
sult for lattices of large temporal extent, T 00, whilst 
here we need to carry out the perturbation theory for 
varying finite T = 3L. We describe the finite volume 
perturbation theory for the NRQCD evolution equation 
in appendix [Xj It turns out that a minor modification 
of the rules for automated Feynman perturbation theory 
account for the effects of finite T in the one-loop case. 

For t large enough, we have that 



G{p,t,(},L) = ^^e(^°+P'/2*^-'=+-)*, 



(46) 



and by fitting to this form for a range of values of p we 
can, in principle, extract the renormalization constants 
Z^f), Z^Q and Eq. However, for the current work we do not 
need Zm^ as we extract Mpoie using Eq. rather than 



Eq. (§ since, as remarked in section [ll Ap the statistics 
available are not sufficient to extract a reliable value for 
Zmo- We therefore evaluate G for p = and measure 
Eq{(3, L), the energy as a function of /? and L. 

From the boundary condition we have G(p, t — 
0, P,L) = 1 and so we cannot fit to the asymptotic form 



below some value t = t^, 



It is a feature of Coulomb 



gauge that Z^ is very close to unity. This is borne out by 
our one loop perturbation theory and also by simulation. 
Consequently, considering Z^ and Eq as functions of t, 
we expect the t dependence of Z^ to be small compared 
with that of Eq and that tmin is not too large . Whilst 
accounting for the need to measure in the asymptotic re- 
gion by fitting only for t > tmin it is useful to account 
for any residual t dependence by including a transient 
function of t in the exponent in Eq. ( 46 ) . From the finite 



volume perturbation theory and from Eqs. (All I and 



(|A12|) E^o^\L,T,t) and Z^^\L,T,t) depend on t, and a 



fit to' their t dependence for small t gives a good indica- 
tion of the explicit transient function we should choose. 
Using the one-loop calculation in this way, we find that 
to extract Eq(P, L) from the high-/? simulation the form 
for G should be chosen as 

G{Q,t,p,L) = Zv,(/?,i)e(^''('''^)*+'^/'), t > Umn. (47) 



G(p, /?, L) = ( V Re Tr {Tle-'P ''G{x, p, t)) 



LJi 

(45) 



where, in practice, we choose tmin — 5 for all L. 

We fit Eq{/3pi,L) to a joint power series in ay^\q*) 
and 1/L, with Uf = 3. In order to do this we need to 
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compute the value of ay (q*) given the value of /Spi with 
which the quenched configurations were generated. We 
first compute ay^ {q* ) from the measured plaquette using 
perturbation theory. The lattice coupling a^, deduced 
directly from the value of /Spi , can be expressed as a per- 
turbation series in a^-'\q*) for any Uf. We eliminate 
by equating the series for nf ~ with that for Uf — 3 

and so deduce a power series for ay\q*) expanded in 

powers of ay\q*). In this way we compute the required 

value of ay (q*) for each value of /3p/. The details follow. 

We choose the ^-scheme defined in terms of the colour 
Coulomb potential and the value of q* is found by us- 
ing the BLM procedure [221 IH] applied to the heavy 
quark self-energy for determining Eq; Miiller 06] gives 

q* = 0.794a~^ for this case. To determine ay\q*) given 
/3 we use the value of the Wilson plaquette, 

from our configurations to calculate cty •* {q* ) using the 
perturbative expansion of Wn. The BLM procedure 
gives the optimal value of q* = 3.33a~^ for this quan- 
tity [1011171 HH]. Note that we compute a^y\q*) in this 
manner, i.e. for Uf =0, since we are using quenched 
configurations. Then we have (nj — 0) 

log(T^ii) ^ -3.068a^^(g*) - 0.5945(2)a^'(q*) 

-0.589(38)a[!?^(g*)2 + . (48) 

We do not find any dependence of Wn on L since it is a 
short-distance, UV, quantity. We now relate ay\q*) to 
aL{a) using [ffil|49] 

ada) = a^;;'\q)(l-vt'\q)a^^^\q) 

-v'^'\q)a'y''\q?) . 

v^^^\q) 2^0 log(7r/g) -1-3.57123- 0.001196n/, 

v^^^\q) = 2/3ilog(7r/g)- [w|"^^]2-f 5.382 - l.OSlln/, 

where /3o and /3i are the coefficients in the /3-function: 
12 1 38 

and then use this expansion to re-express the result in 
terms of ay'\q*). We find 

ap\q) = af{q){l + m{q)a^^\q)+u^{q)a^^\qf), 

Mq) = v[''^\q)-V^f\q), 

U2{q) = V^^''\q)~vi'\q)+u,iq)v^^'^\q). 

We then run av{q*) from q* = 3.33a^^ to q* = 0.794a~^, 
appropriate for the fit to Eq(P,L), using the three-loop 
running 

dbg/i'^"^ = -av{^J-)^ {l3o + f^iavifJ-) + |32vav{^J'f) , 
P2V = (4224.18 - 746.006nf + 20.8719n?), 



where we suppress the n / superscript from now on, using 
rij = 3 implicitly. 

We fit G{0,t, (3, L) separately, as discussed above, for 
the set of /3, L values and deduce Eo{l3,L). As the 
data may contain residual auto-correlations, we resam- 
ple via blocking to determine the true statistical error. 
Within independent chains, sequential measurements are 
grouped together into bins and the means of each bin are 
treated as statistically independent. The size of the bins 
is determined by examining the scaling of the variance as 
a function of the bin size, and is dependant on the values 
of L and /3pj, and the operator being measured. We then 
fit these values to the form 

EoiP,L) - iE^o'\L,T/2) + S)av{q*) 

+ (c20 + ^C2i)av{q*)^ + C3oay(g*)^ (50) 

with q* = 0.794a-i and T = 3L. Here E^q^\l,T/2) is 
the calculated value for the one-loop contribution which 
includes the contribution from tadpole improvement of 
the NRQCD Hamiltonian; this contribution is a constant, 
independent of f3 and L. We allow for a small additive 
adjustment S, independent of /3pi and L, in the values of 

the Eq^\l,T/2) accounting for any minor mismatch be- 
tween their analytical and numerical calculation; as we 
should expect, S is found to be very small. The finite-size, 
L, dependence of Eq is included in Eq^\l,T/2) and in 
the two-loop coefficient. We find that this parametriza- 
tion is sufficient for a very good fit to the data; within 
errors we do not discern any o? jl? or o? jL contribu- 
tions. The fit is for 116 degrees of freedom (4 parameters, 
15 /3 values and 8 L values) and we find = 1-2 and 
1.1, respectively, for aniQ — 1.72,2.5. In Fig [2] we show 
EQ{f3, L) plotted versus av{q*) for the different L and for 
am-o = 1.72. The quenched results that we require are 

-C'O — C20- 

IV. NONPERTURBATIVE DETERMINATION 

OF £,i„. 

We now briefly discuss the nonperturbative determina- 
tion of the meson energies iS'sim- The method is standard 
and this NRQCD action [T7] has been thoroughly tested 
by HPQCD in a range of calculations. 

We use two ensembles of gauge configurations gener- 
ated by the MILC collaboration with ny = 2 + 1 ASQtad 
sea quarks, which we denote coarse (^0.12 fm) and fine 
(~0.09 fm) t25i^. Details are given in table |ll) The 
light quark masses on these ensembles are not particu- 
larly chiral but we have seen that the light sea quark 
mass has negligible effect on most quantities in the bot- 
tomonium spectrum [TTj. The lattice spacing on these 
ensembles has been determined using the static quark 
potential parameter ri in |50j . and is given in the table. 

The NRQCD action is given in section |II C| and in- 
cludes one-loop radiative corrections to the coefficients 
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0.08 




FIG. 2: Eo{av(q*),L) for aM = 1.72 for both data and fit 
for the values of lattice size L"^ x T , T = ZL used in the 
extraction of the two- and three-loop quenched coefficients in 
the perturbation series for Eq. We write Eo as a function of 
the av{q*) value rather than jSpi. Here, q* = 0.794a~^. This 
fit has = 1-2. 



TABLE II: Details of the two ASQtad gauge configurations 
used in the nonperturbative determination of iJsim- /3 is the 
gauge couphng, a^^ is the inverse lattice spacing determined 
using the static quark potential parameter ri, uoami, uoarus 
are the light sea quark masses, L and T are the lattice dimen- 
sions and ncfg the size of the ensemble. 



Set 


/3 


a-i (GeV) 


uoami 




LxT 


"•cfg 


coarse 


6.76 


1.652(14) 


0.01 


0.05 


20x64 


1380 


fine 


7.09 


2.330(17) 


0.0062 


0.0310 


28x96 


904 



calculated in [T71 [5T] . The coefficients are listed in table 
III The same coefficients are used in the perturbative 
calculations and in the high-/3 simulations, but with as 
evaluated at a scale appropriate for /?, as discussed in 
section III CI 

Tuning the bare 6-quark mass accurately is an impor- 
tant part of the calculation as this is a potential source 
of error in TOh(mb)- The heavy quark energy shift means 
that we cannot tune using the meson energy directly but 



TABLE IV: Masses and extracted energies from the nonper- 
turbative simulations, amo and arus are the bare (valence) 6 
and s masses, aEsim.x are the fitted ground state energies of 
the meson X in lattice units. The first row is for the coarse 
ensemble and the second for fine. The errors are from statis- 
tics/fitting only. 



amo 










aEsiin,B* 


2.50 
1.72 


0.0496 
0.0337 


0.46591(6) 
0.41385(4) 


0.42579(3) 
0.38124(2) 


0.6278(5) 
0.4812(5) 


0.6595(6) 
0.5027(7) 



we must use the kinetic mass determined from the dis- 
persion relation, which is much noisier. A detailed study 
of the systematic errors incurred and their effect on the 
accuracy of the bare mass was carried out in Ref [T7] . To 
reduce systematic errors we use the spin average of the 
vector and pseudoscalar bottomonium states 



^ (3Mki„,T + A/ki„^„J/4, 



(51) 



which eliminates errors from missing spin dependent 
higher order terms and radiative corrections in the ac- 
tion. We must also take account of missing electromag- 
netic effects, sea charm quarks and annihilation of the 
Tjb to gluons by shifting the experimental values appro- 
priately. These effects were estimated in |20|, resulting 
in an adjusted experimental value of M^^"^^ — 9.450(4) 
GeV, where the error comes from taking a large un- 
certainty on the shifts that were applied. The cor- 
rectly tuned bare 6-quark masses in lattice units that 
we obtain are 2.49(2)stat(l)sys on the coarse lattice, and 



1.71(2),tat(l)s 



The first error includes a sizeable sta- 



tistical error from the kinetic mass and all lattice spacing 
errors, the second includes the systematic errors in the 
kinetic mass estimated in [T7]. The effect of these errors 
are included in the final error budget. 

The valence strange quark propagators used in the 
Bg mesons use the Highly Improved Staggered Quark 
(HISQ) action |52] and are tuned using the rjs meson. 
This is a fictitious ss particle which, with the addition 
of experimental data for M^r , Mk and chiral perturba- 
tion theory, is a very convenient choice for tuning the s 
mass and fixing the scale. The value on the = 2 -I- 1 
ensembles that we are using is M,,^ = 0.6858(40) GeV 



TABLE III: Coefficients used in the nonperturbative simu- 
lation. uo,p is the plaquette tadpole improvement factor, Ci 
are the coefficients in SH. 



Set 




Cl 


C2 


C3 


Ci 


cs 


C6 


coarse 


0.86879 


1.31 


1.0 


1.0 


1.2 


1.16 


1.31 


fine 


0.878214 


1.21 


1.0 


1.0 


1.16 


1.12 


1.21 



The ground state energies E'sim are extracted from 
multiexponential Bayesian fits |53| to meson correla- 
tion functions that use multiple smeared sources for the 
quark propagators. To further improve statistics we used 
stochastic noise sources and ran 16 time sources on each 
configuration for the T, and 4 for the Bg. The results 
are listed in table lIVl 
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V. CALCULATING THE MS &-QUARK MASS 

Now that Eq and -Esim have been determined, we can 
combine the results into a perturbative series for fnf){fnb) 
in the MS scheme. This requires various scheme conver- 
sions and changes of scale to give the series at the scale 
relevant for the fo-quark mass. This then gives the result 
at n/ = 3 and we can use known formulas to convert this 
to the usual Uf = 5 result. We repeat this whole process 
at both values of the bare mass to check for discretisation 
errors which will then be included in our error. 

To further reduce systematic errors, we adjust Eq [T] 
so that we use the spin-averaged bottomonium mass 
Mjji = (3Mx + -^?7b)/4. This removes any error from 
spin dependent terms in the NRQCD action. As dis- 
cussed in section [TV] the experimental result used must 
be adjusted to M^6,oxpmt = 9.450(4) GeV to reflect the 
absence of electromagnetism, sea charm quarks and rjh 
annihilation. 



A. Perturbative series for mb{mb) 

So far, all our perturbative results have been expressed 
in terms of ay, the coupling constant defined in the V- 
scheme at the scale q* = 0.794/a. 



uEq — aEQ^^aviq* 



aE, 



(2) 



aE, 



uOJ 



+ aE^'^''^aUq* 



(52) 



The results for each component are given in table |V] 

The series expansion of uEq is truncated at and we 
take Uf — S as this is the number of sea quarks in the 
nonperturbative determination of Egim- No fermionic 
contributions are included in the series. The effects of 
the one loop tadpole corrections are directly included in 
the tadpole improved results from the high-/? simulation, 
as are the quenched two loop tadpoles. However, the 
two loop fermionic tadpole contributions are not included 
in the high-/? results so we must add the corresponding 
correction, gEq'^'-^ , to the energy shift. aE^'^'^ is given 
by [3H] 



aE, 



uO,/ 



2amo 



,(2),/ 



(53) 



(2) f ' 

where Uq is the fermionic contribution to wo,p given 
in section IIII Al 

The other perturbative factor that we need is the pole 
to MS renormalization Zm which is reproduced in ap- 
pendix [C] Inserting these two series into eq. ([t]) gives a 
series for mb(TOb). 

We now relate av{q*) to aj-ig{q*). This is done using 
the three-loop relation in [54M56] which is summarised 
in appendix O and express Eq as a series in the MS 
scheme. Matching is done at q* to avoid logarithmic 



contributions. The series is then run to /i = 4.2 GeV 
using the 4-loop MS beta function. 

To evaluate the series we need the relevant value of 
ajYg, which in this case is the 3-flavour value at mfc. 
Since MS is a mass independent scheme, high mass 
particles do not explicitly decouple from the beta func- 
tion and one must construct an effective theory with 
ni — nj ^ 1 quarks when crossing a quark mass thresh- 
old [57J. This introduces discontinuities in the running 
of ajjg at the thresholds which have been calculated 
to 4- loops in [SH], and we give the relevant formulas 
in appendix [C] We start with the current PDG average 
ajYs{Mz,nf = 5) = 0.1184(7) which we run to 4.2 GeV 
using the 4-loop running with nj = b |59j , then matching 
to the Uf — 4 theory and running down to 1.2 GeV to 
match to Uf = 3, before running back up to 4.2 GeV with 
Uf — 3 running. We find ajjgirfii^nf = 3) = 0.2159(20). 
Small changes in the matching scales have negligible ef- 
fect on the value. 

Using this value of the coupling the results using Mi^^ 
are fnii{rni,,nf = 3) = 4.195(8) GeV on the coarse lattice 
and rni,{rni,,nf = 3) = 4.198(10) GeV on the fine lat- 
tice. We also tried allowing the scale to float and solving 
such that fi was exactly the MS mass but this makes 
negligible difference to the result. The results using the 
Bs mass give fni,{fnb,nf = 3) = 4.177(8) GeV on the 
coarse lattice and rnb{fnb,nf = 3) = 4.191(10) GeV on 
the fine lattice. These are consistent with the bottomo- 
nium results. This error includes statistical errors in the 
perturbation theory integrals, lattice spacing error, and 
simulation errors in the ground state masses (negligible) . 
We have not yet included an estimate of the truncation 
error in the perturbative series. 

Our calculations were performed using lattice results 
with Hf — 3 sea quarks. In order to compare to the real 
world we must match this value to n^^ = 5. As with 
the coupling constant, a running quark mass in a mass 
independent scheme is discontinuous at flavour thresh- 
olds and must be matched to an effective theory with a 
different number of flavours. The formula for the mass 
decoupling is given in appendix [C| in equation ( C12 ). We 
run down to 1.2 GeV with three flavour mass running 
[59], match to a theory with nf ~ A, run up to 4.2 
GeV and match to the Uf = 5 theory. Again, small 
changes to the matching scale or the final scale at which 
we evaluate the mass have negligible effect. After this 
running, the values we obtain for the M^b results are 
rnb{ffib,nf — 5) = 4.161(10) GeV on the coarse lattice 
and nib {rrTb,nf = 5) = 4.164(12) GeV on the fine lat- 
tice where from now on we state rt/ explicitly. Overall, 
matching to the ny = 5 theory shifts the mass down by 
around 30 MeV. 

In principle there may be discretisation errors arising 
from lattice artefacts. Since we have two lattice spac- 
ings available we can fit the results as a function of a 
to obtain the physical result and to allow a systematic 
error for this dependence. In fact the dependence is very 
mild as is clear from the fact that all of the results are 
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TABLE V: Perturbative results required to extract the MS 
mass. The quenched results, indicated by superscript q, are 
from high-/? simulations. The one-loop data are the exact 
pertiu'bative results extrapolated to infinite lattice size. The 
two loop results include both quenched and fermionic contri- 
butions. The three-loop values include only quenched results. 
We evaluate all results in the V-scheme at a characteristic 
scale of q* = 0.7940"^ [46] . 



amo 



aE, 



(1) 



aEl 



(2) 



aE, 



uQ,f 



aE, 



(3),9 



2.50 0.6786(1) 1.16(4) 
1.72 0.5752(1) 1.30(4) 



0.2823(6)n/ 0.158531(16)n/ 2.3(3) 
0.3041(3)n/ 0.186607(19)n/ 2.3(3) 



consistent with each other. Our NRQCD action contains 
discretisation corrections that get renormalized as a func- 
tion of the cutoff amo £^nd so we allow an additional mild 
dependence of the fit function on amo. This makes no 
difference to the fit. The form is 



mb{mh){a,Sxm) = mb{mb) 



1 + ^ dj(Aa)^-' (l + djbSx„i + djbb{5xmY) 



(54) 



where we have allowed discretisation effects with a scale 
of A = 0.5 GeV and cutoff dependence via 5xm = (amo — 
2.1)7(2.5 — 1.7) which varies between ±0.5. Priors on the 
values are 4.2(5) for the mass, 0.0(3) for the a? term since 
our action is one-loop improved, and 0(1) for everything 
else. 

Some of the errors in the data are correlated and we 
allow for this in the fit. We multiply the Wib values by 
a (1 -|- JT-z^s) truncation error (discussed below) which 
is 100% correlated between the points on the two lat- 
tice spacings. The errors on all quantities coming from 
the high-/? simulations are correlated with correspond- 
ing errors on the other lattice spacing. Statistical errors 
coming from VEGAS integrals are uncorrelated. 

We only fit the bottomonium results as the Bs re- 
sults are in very good agreement. The result of the fit is 
nib {mb,nf = 5) = 4.166(42) GeV. 



B. Error budget 

Broadly, the three main sources of uncertainty in our 
result for the 6-quark mass are: statistical errors, errors 
from truncating the perturbation series and other sys- 
tematic errors. We expect the 0{al) perturbative con- 
tributions to dominate the uncertainty in our final result. 
In this section we discuss each of these sources of error 
in turn and tabulate our error budget in Table |VI[ 

a. Statistical errors Statistical errors arise in the 
nonperturbative calculation of -Bsim , and in the contribu- 
tions at each order in the expansion of the heavy quark 
energy shift, Eq. The statistical error in Egim comes from 
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FIG. 3; Results for the n/ = 5 MS mass using both bot- 
tomonium and Bs meson simulation data, and the fit to the 
bottomonium results. The errors on the data points include 
statistics, error on Qjyjg and a correlated truncation error on 
the perturbative series. Additional (subdominant) errors are 
described in the text. 



the fit to lattice 2-point functions and is completely neg- 
ligible. The statistical error in the one-loop piece of Eo 
comes from the evaluation of diagrams using VEGAS and 
from the extrapolation to infinite volume. The uncertain- 
ties in the two loop and three-loop quenched coefficients 
of Eg arise from the simultaneous fit to a and L. This 
is significant at 14 MeV. The statistical error in the two 
loop fermionic coefficient is due to the numerical evalua- 
tion of the Feynman diagrams and the extrapolation to 
zero light quark mass. 

b. Perturbative errors The three-loop fermionic 
contribution to the energy shift is unknown, so we esti- 
mate the error due to this contribution as 0{nf x ctj^)- 
This is the dominant source of error in our calculation. 
Perturbative errors from running the coupling and quark 
mass are negligible as the formulas are higher order. 

The fermionic contributions are the only unknown 
source of uncertainty at three- loops in our result. In 
principle these effects can be calculated using automated 
lattice perturbation theory. However, there are a large 
number of diagrams to evaluate, many of which are likely 
to have complicated pole structures and possible diver- 
gences (the energy shift is infrared finite, but individual 
diagrams may have divergences that ultimately cancel). 
The complexity of such a calculation would be consider- 
able. 

c. Other systematic errors 

• Bare mass tuning: The tuning of the bare 6-quark 
mass used in Eq and i?sim is a source of error. We 
can estimate the error due to mistuning using the 
errors given on the tuned masses 2.49(2)stat(l)sys, 
and 1.71(2)stat(l)sys and by estimating the bare 
mass dependence of each quantity. We use only 
the one-loop piece of Eq and compute the value at 
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an extra mass, we find a linear dependence with 
a slope of 0.13. For i?sim, we use the results at 
different bare masses given in [T7] and find a de- 
pendence that is less than 0.01, which we take to 
be linear for these small increments. By recomput- 
ing fnf){rnb) taking a la deviation in the bare mass, 
we find errors of 4 MeV on the coarse lattice, and 
6 MeV on fine. We take the larger of these as an 
error on our result. 

• Corrections for missing electromagnetism, charm 
quarks in the sea and r/b annihilation were es- 
timated and applied to the experimental T,rib 
masses. We add the errors linearly rather than in 
quadrature and propagate this error through to the 
final result, which gives 1.9 MeV. 

• Higher order relativistic corrections: These arise 
from not including 0{v^) terms in our NRQCD ac- 
tion and, with ~ 0.1, could contribute 1% of the 
binding energy which is 5 MeV. 

• Radiative corrections: a'^v^ should be smaller at 
around half a percent of the binding energy so we 
take 2.5 MeV. 

• Lattice spacing errors, including ri/a: These are 
included as the "statistical" error on the data 
points in the plot but we estimate their contribu- 
tion to the final error to be 4.5 MeV. 

• Lattice spacing dependence: We incur an error 
from fitting the two masses as a function of a which 
we can estimate from the fit. The lattice spacing 
dependence is not significant but we find 16 MeV, 
this is already included in the total error quoted 
from the fit. 

• Sea quark mass dependence: We have only used 
one sea quark mass in our calculation but in pre- 
vious calculations we have observed very mild de- 
pendence in Esiai |17j. Errors from light sea quark 
mass dependence should be negligible compared to 
our other errors. 

With these errors included, our final result for the MS 
&-quark mass is: 

mb{mb,nf = 5) = 4.166(43) GeV. (55) 



VI. DISCUSSION 

We can compare our result to previous values from the 
literature. As discussed in Section |T] there are a num- 
ber of accurate theory results from comparing contin- 
uum QCD perturbation (through a^) for moments of the 
vector charmonium current-current correlator to exper- 
imental results extracted from a{e'^e~ — >■ hadrons) in 
the b region. In [T], for example, the result m,b{mb) — 



TABLE VI: The 6-quark mass error budget, systematic error 
estimates are discussed in more detail in the text. 



Source 


Jirror (iVieVJ 


Jirror (/o) 


UfCxi, perturbative error 


36 


0.9 


Mr , Mrii, experiment 


< 0.1 


< 0.01 




< 0.1 


< 0.01 


amo tuning 


6 


0.14 


VEGAS integration 


< 0.1 


< 0.01 


High-/3 statistics 


14 


0.35 


a dependence 


16 


0.38 


Scale uncertainty 


4.4 


0.10 


Os uncertainty 


0.2 


0.01 


Relativistic 


5 


0.12 


Radiative asU* 


2.5 


0.06 


E&M, Charm sea, annih. 


1.9 


0.05 


Total 


43 MeV 


1.0 % 


4.163(16) GeV is obtained. 


In [S] lattice 


QCD calcu- 



lations of time-moments of the r]b correlator are used 
instead of the experimental results to give mbirfib) = 
4.164(23) GeV. It was important in this calculation to 
use pseudoscalar correlators in a lattice QCD formalism 
(HISQ) that has absolutely normalised pseudoscalar cur- 
rents. Our result agrees with these two values. It is not 
as accurate because we are not using such high order 
QCD perturbation theory but it nevertheless provides a 
check from a completely different perspective at the level 
of 1%. 

There are also a number of results using alternative 
methods from lattice QCD but these are not typically 
very accurate. An early result for nib with NRQCD b- 
quarks on the n/ = 2 + 1 MILC configurations including 
M, d and s sea quarks was 4.4(3) GeV [8^, the large error 
here arising from the use of one-loop lattice QCD pertur- 
bation theory for Zm ■ More recently, methods have been 
developed by the ALPHA collaboration for determining 
the energy shift for lattice Heavy Quark Effective Theory 
nonperturbatively, including next-to-leading-order terms 
in the inverse heavy quark mass expansion for the va- 
lence &-quarks |60| . This has been implemented on gluon 
field configurations including u and d sea quarks in the 
clover formalism. Combining with the experimental B 
meson mass in a similar approach to the one used here, 
gives mbirfib) = 4.22(11) GeV. The error here is dom- 
inated by lattice statistical and systematic errors. An- 
other method by the ETM collaboration (il uses a ra- 
tio of quark masses to heavy-light meson masses with a 
known infinite mass limit. This is implemented on gluon 
field configurations including u and d sea quarks in the 
twisted mass formalism and valence b and light twisted 
mass quarks. Interpolating to the 6-quark and using ex- 
perimental meson masses gives: Wib{rrib) = 4.29(14) GeV, 
with an error dominated by lattice statistical errors. Note 
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FIG. 4; Comparison of our result with other recent theory- 
based 6-quark mass determinations. We include all determi- 
nations listed in the PDG summary table [62] but separate 
lattice QCD determinations with nj = 2 and nj = 2-1-1 sea 
quarks for easier comparison [T] [21 151 [SHI [HIl . 



that neither of the ALPHA or ETM results include s 
quarks in the sea and the error from this is not estimated. 

Fig. |4] collects a number of lattice and continuum QCD 
determinations of the 6-quark mass for comparison. The 
evaluation of 4.18(3) GeV in the Particle Data Tables [51] 
is shown by the grey band. There is good consistency 
between all determinations including the new result of 
this paper. 



VII. 



CONCLUSION 



In this paper we have presented a new determination 
of the 6-quark mass from simulations of lattice NRQCD 
at two heavy quark masses. The uncertainty associated 
with previous determinations of the &-quark mass from 
lattice NRQCD was dominated by the one loop pertur- 
bative calculations used to extract the 5-quark mass. By 
calculating the heavy quark energy shift at two loops, 
we have significantly reduced this uncertainty. The re- 
sulting error is now in line with the most precise lattice 
determinations available. 

In order to efficiently calculate renormalisation param- 
eters at two loops, we implemented a mixed approach, 
combining quenched high-/? simulation with automated 
lattice perturbation theory. We were also able to ex- 
tract estimates of the three loop quenched contributions 
to the energy shift from high-/? simulations and found 



that all perturbative coefficients are well-behaved. The 
reliable extraction of the two loop energy shift convinc- 
ingly demonstrates the effectiveness of our approach. 

As part of this calculation, we also determined the 
fermionic contributions to the two loop tadpole improve- 
ment factor for both the Landau and plaquette tadpole 
definitions. 

We undertook a number of checks of both the auto- 
mated lattice perturbation theory and the high-/? simu- 
lations. For the former, we confirmed that we could re- 
produce published one loop results, that the energy shift 
was infrared finite and that the fermionic insertions in the 
gluon propagator obeyed the relevant Ward identity. For 
the latter, we were able to compare one loop results to 
the exact finite size perturbation theory results to ensure 
the correctness of our fits. 

The uncertainty in our result is now dominated by the 
unknown fermionic contributions to the three loop energy 
shift, which is in principle calculable with automated lat- 
tice perturbation theory. Greater statistics in the high-/? 
simulations may also allow us to extract the quenched 
contributions to the mass renormalisation with sufficient 
precision to enable an independent determination of the 
6-quark mass by direct matching. 
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Appendix A: Finite volume perturbation theory 

Without loss of generality we consider a scalar model 
that is sufficient to demonstrate the approach. We take 
the NRQCD evolution for the heavy quark Green func- 
tion to be 

G(p,i) = i^(p,t-l)G(p,t-l), (Al) 



where 



K{p,t) = ifo(p,i)(l-|0), 

/ 2 \ " 

K,{P) = 1-;^^ • (A2) 
\ Iran J 
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Here K{x, t) is the approximation to the evolution oper- 
ator with 



H=\ 1 



V2 

2m 



We then have that 



GQ{p,t) = K{p)\ with Go(p,0) = l. 



(A3) 



(A4) 



The diagram we consider is the rainbow diagram. The 
vertices are labehed with (p. t) coordinates, appropriate 
to the Hamiltonian formahsm. The vertices are separated 
by time r. The rainbow diagram has r > whilst the as- 
sociated tadpole diagram has t — Q. There is no effect of 
finite T on the calculation of the tadpole diagram, which 
is therefore given by finite-L Feynman perturbation the- 
ory. At 0{g^) from the diagram we have the contribution 

T-l 

G2{p,t)=g^ E {t-T)Ko{pY-'T{q,T)Ko{p-qr. 

q,r— 1 

,(A5) 

The factor {t — r) is the number of temporal positions 
the graph can adopt and T{q, r) is the (^-field propagator, 
given by 



T-l 



r(9,r) = ;^^r(q,9o)e^' 

r(q,go) 



S1=0 



1 



(A6) 



where 



90 



9t = 



2 ttQ, 
L ' 



(To = 2 sin 



go 



(A7) 



qi = 2 sin 



with <V, <T and < Qi < L. Then the contribution 
from the rainbow diagram is 



t„2 



G2ip.t)=Koipyg 



xr(g,(7o) 



nr=i 



We now let 



R{p,q) = 



Koip) 

Ko{p ~ q) 
Ko{p) 



(A8) 



(A9) 



where q = (qo^q). Then, using Eq. ( |A4[ ), the one-loop 
rainbow diagram correction to the Green function is 

G{p,t) - [l-g^A{p,t)t + g^B{p,t)]Goip,t) 

~ [l + g2B(p,t)]e-9'^(p.*)*Go(p,i), (AlO) 



and we deduce that 



g'^A{0,t) and 

zl^\L,T,t) = g^B{0,t) 
depend on t but that for t sufficiently large both quanti- 
ties will approach their asymptotic value. We then have 



zl^'\L,T,t) ^ g' 



, 1 



^^r(q,go)i?(0,'Z)^(All) 

Qi:^ T = l 
t 

5]^rf(q,go)i?(0,9r-(Al2) 

Qi,n T=i 



We first consider Eq^\l, T, t). We carry out the geomet- 
rical sum and find 



A{p,t) = -g^^y f(g,go) . 



(l-i?(p,g)*). 



For |i?(0,g)| < 1, for q, the limits T 
be taken. We have that 



oo, t 



(A13) 
oo can 



R{p,q) 

l-R{p,q) 



Ko{p-q) 
e^iPo-qa) - Kq{p - q) 

e»(Po-9o)^^(p _ g)Go(p - q),(A14) 



where we have used the on-shell condition for the external 
quark: e*^° — Ko{p). In this limit we find 



A{p, oo) 



1 



-g 



E 



dz 



-i{po-qo) 



\z\ = l 



xKoip^q)GQ{LU,p-q) , (A15) 

where w = e*^^''^'^''^ with z = e^*'°, and the integral 
is over the unit circle in the complex z-plane. This is 
the expression for the rainbow diagram derived from the 
NRQCD Feynman rules applicable in the limits T ^• 
oo, t — oo. 

We conclude that to account for the effect of finite 
temporal extent of the lattice in the perturbation theory 
we make the replacement 

Go(c^,p - q) Go(c^,p - q)[l - R{p, qf] (A16) 

for the internal quark propagator, and carry out the sums 
over the discrete values of q and qo- There remains the 
choice for the value of t in this expression. We found 
that the results were insensitive to this choice as long as 
t was not close to either or T and so we chose t = T/2 
for our calculations. R{p, q) is computed automatically 
by a numerical search for the poles of the external and 
internal propagators which gives Kq (p) and Kq [p — q). 
For e'^^'^ (L, T, t) we set p = 0. 

The wavefunction renormalisation, .^i^'', is given by 



Z^'\L,T,t)^^g' 



——yy 



r(<7,go) 



(A17) 



15 



evaluated on-shell: e"*''" = ifo(p). This is the usual 
formula applied to our augmented Feynman rule and 
the derivative is computed using our automated taylor 
derivative procedure. 

In some cases we can have \R\ > 1. This is the situ- 
ation for some values of q given p and certainly occurs 
in moving NRQCD (mNRQCD) [75]. Because NRQCD 
is in the Hamiltonian formalism the value of t in Eq. 



( A8 ) is finite and the singularity in the quark propagator 
is removable. The poles in the gluon propagator are at 
z = z± with \z±\ ^ 1 and z+Z- = 1. Schematically, Eq. 



(A15) takes the form 



dz- 



1 z(l-(a/z)*) 



(A18) 

where C is a constant and a — Ko{p — q)/KQ{p). The 
integration contour is \z\ — 1 and is determined by the 
formalism; no distortion is available in the NRQCD evo- 
lution to avoid pole crossing. However, the singularity at 
z = a is removable and so there is no issue of it cross- 
ing the contour. The integration is done by Cauchy's 
theorem at the z = z^ pole, and the factor from the geo- 
metric summation is then evaluated to be (1 — (a/z+)*); 
the need to consider the pole of order {t — 1) at the origin 
is then avoided. Since \a\ < \z^\ the limit t — >■ cxd can now 
be taken. This corresponds to the usual rule for analytic 
continuation in the calculation of the Feynman diagram 
where the radius |z| of the contour is increased to avoid 
crossing by the quark pole at z = a. 



Appendix B: Generating Configurations and Gauge 
Fixing 

1. Langevin Marltov Chain Configurations 

Configurations for the Monte Carlo simulations are 
generated with a Markov chain that is updated via a 
Langevin algorithm. The Langevin method treats the 
Markov chain as a classical path in phase space, using 
the action as a potential to enforce the Boltzmann distri- 



bution. Using the notation of section IIIB the Langevin 
equation is given by 



dU 



dS 
dU 



(Bl) 



where S is the action, 77 is a random noise term, and r 
is the distance along the path. Using the Fokker-Plank 
equation it can be shown that this path will sample the 
configuration space with probability density 



PiU) 



-siu] 



(B2) 



the Boltzmann distribution, as desired. 

As Eq. (Bl) is an initial value problem, its solution 



can be approximated via an iterative method, where the 



derivative on the left hand side is written as a finite dif- 
ference, with step size e. This introduces step-size errors 
in the action so that the distribution that is simulated is 
altered to 



P{U) = e" 



-S\U,e] 



(B3) 



where S[U, e] is the simulated action which is expansible 
as 



S[U, e]^S + eSi 



e^S2 



(B4) 



The step-size errors in Eq. (B4) can be systemati 



cally eliminated using higher order approximations to 
the derivative in Eq. (Bl). In this work a second or- 



der Runge-Kutta algorithm (RK2) which eliminates 0{e) 
errors was used. This is implemented as a mid-point 
method adapted to diffusion on a group manifold |77j . 

Simulations were run with a step size e = 0.2. Analy- 
sis shows that e scaling errors are of the order w 0.05%. 
Auto-correlation times were measured to be of the or- 
der of 5 — 10 (25 — 50 updates), for the plaquette and 
10 — 20 (50 — 100 updates) for the twisted Polyakov loop 
[To] . Here 100 configurations were skipped between mea- 
surements. For each value of /3pi on each lattice size, 32 
independent Markov chains were generated. Each chain 
produced 128 configurations (4096 configuration in to- 
tal). 



2. Gauge Fixing with Twisted Boundaries 

Configurations generated from the Markov chain have 
the gauge freedom described in Eq. (24 1. This can be 
fixed by the application of a gauge condition. In this 



work we wish to fix the configurations to Coulomb gauge. 
In the continuum Coulomb gauge is achieved by the the 
gauge transformation that satisfies 



a.Af = 0. 



(B5) 



On the lattice this corresponds to maximizing the quan- 
tity 



x,i—l 

1 

Te 



g{x) Ui{x) Ui{x + Bi) g'^{x + 2ei) 



, (B6) 



with respect to the gauge transform field g{x) for each 
time slice. This is 0{a^) improved [7S]. The maximisa- 
tion is preformed via a conjugate-gradient method, using 
a backtrack line search. Each time slice is gauge fixed 
separately. Errors due to numerical maximisation are 
estimated to be insignificant. 

Fixing to Coulomb gauge leaves an ambiguity, since 
it is possible to construct an additional purely temporal 
gauge transformation 

Ui{x,t) ^ Uf \x) = g(^)(i) Uiix)g^^^Ht + l). (B7) 
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This gauge transformation must obey the twisted bound- 
ary conditions 



for i = 1,2,3. The only solutions are 



(B8) 
(B9) 



for n — 0, . . . , Nc, where z is given in ( p2| ). 

After fixing to Coulomb gauge, each time may be in 
a different gauge. In order to measure time dependent 
operators, the time slices must all be in the same gauge. 
The gauges are all fixed to be the same as that on the 
first time slice. Since the gauge transformation in (B9) 



form a group, this is achieved by applying an additional 
transformation. The gauge transformation on the first 
time slice is set to the the unit matrix 



1. 



(BIO) 



The transformations on subsequent time slices are chosen 
sequentially to maximise 



Re Trg'^^\t-l)U4{0,t-l)g 



for i = 1,...T- 1. 



(T)t 



(Bll) 



Appendix C: MS matching formulas 

The relation between ay and aj^g is given by 



The pole to MS renormalization is calculated to three 
loops in [16] 



M, 



pole 



with 



ZM{m.b)mb{mb), 



(C6) 



7 (rn\-^ '^Msi^b) 



3 IT 
2 



(-1.0414n/ + 13.4434) 



(0.6527n^ - 26.655n/ + 190.595). 

(C7) 



We actually need the inverse of this series which we define 
as the 3-loop approximation to 1/Zm- With = 3 this 
is 



Zui^b) = 1 



G.42441318a 
0.86542701a- 



MS 

2 

US 



2.94639a|y;g. (C8) 



The MS coupling is discontinuous at quark mass 
thresholds since the heavy mass quarks are explicitly de- 
coupled by matching to a theory with a different number 
of flavours. The formula for matching the Uf theory to a 
theory with n; = n/ — 1 flavours at the threshold is [5S] 



MS 



MS 



f^v = Oi-^{l.Q + coa^ + cia\^). (CI) 
The coefficients are: 



Co = (fli + ^olog(a;))/47r 

ci = (a2 + (/3o log{x)f + (/3i + 2/3oai) log(a;))/(47r)2 

with log(a;) — since both coupling are evaluated at the 
same scale 



/3o 
Pi 

«! 
02 



11 - 2n//3., 
2(51 - 19n//3), 
(31Ca-20T^n/)/9, 
/4343 



(C2) 
(C3) 
(C4) 



V 162 
/1798 



7rV4 + 22C(3)/3 C2 (C5) 



V 81 
55 

y ~ 



-)-56C(3)/3j CaTfUf 
16C(3)) CfTfUf + ^T]n 



'2„2 
/■ 



Note the discrepancy between [3S] and 



with everything evaluated at the threshold scale of the 
rif theory and the coefficients 



C2 
C3 



11 

72' 
82043 

27648 



C(3) 



564731 2633 



124416 31104 



ni. 



(CIO) 
(Cll) 



Crossing thresholds for a running mass in a mass in- 
dependent scheme gives the same difficulties as the cou- 
pling. The relation between the ni flavour effective the- 
ory and the rif flavour theory for the MS running mass 
at the threshold is [TSj 



m 



("0 



m 



(»/) 



0.2060 

^2 



(1.8476 -h0.0247n;) / („^)y 



(C12) 



For the inverse of these operations we include higher or- 
der terms so that it reproduces the original value to bet- 
ter accuracy. 
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